#include<iostream>
using namespace std;
#include<cmath>

/*
 * Nevill 算法
 *
 * 给定 n 个数据点　(x_data[], y_data[])，指定内插/外推点位置　x，Nevill()函数自动进行多项式拟合，并返回 P(x)值
 *
 * 例：
 * P_0
 * P_1	P_{01}
 * P_2	P_{12}	P_{012}
 *
 * 程序定义一个 3x3 的矩阵(虽然上面是一个下三角阵)，然后逐列计算矩阵元的值, 最后返回 P_{012} 的值。
 *
 * 如果数据点的　x　坐标太靠近，误差可能会比较大，所以程序中断。
 *
 */
double Nevill(int n, double *x_data, double *y_data, double x){
	double y;
	//检查数据点　x　坐标，如果太靠近就终止程序
	for(int i=0;i<n;i++){
		for(int j=i+1;j<n;j++){
			if( fabs(x_data[i]-x_data[j]) < 1E-6 ){
				cout<<" datapoints' x values are too close.\n";
				exit(1);
			}
		}
	}

	double ** tempP = new double * [n];//开辟动态内存
	for(int i=0;i<n;i++) tempP[i] = new double [n];

	for(int i=0;i<n;i++) tempP[i][0] = y_data[i];// 第一列：P_0, P_1, \cdots, P_{n-1}
	//use the recurrence relations to derive the final polynomial
	// 用递推公式，推出第二列，第三列。。。
	for(int i=1;i<n;i++){
		for(int j=i;j<n;j++){
			tempP[j][i] = (x - x_data[j])*tempP[j-1][i-1] + (x_data[j-i] - x)*tempP[j][i-1];
			tempP[j][i] /= x_data[j-i] - x_data[j];
		}
	}

	y = tempP[n-1][n-1];

	for(int i=0;i<n;i++) delete [] tempP[i];
	delete [] tempP;//删除动态内存

	return y;
}

int main(){

	int n;//数据个数
	double x,y,dy;

	double *x_data;
	double *y_data; 

	FILE *fp=fopen("Nevill_input.txt","r");
	if(fp==NULL){
		cout<<"error: failed to open Nevill_input.txt.\n";
		exit(1);
	}

	//读一遍输入文件，并将数据行数输出到屏幕上
	n=0;
	while( fscanf(fp, "%*c%*[^\n]") != EOF ){
		n++;
	}
	cout<<"n="<<n<<endl;
	
	//开辟动态内存，用于存放数据
	x_data = new double [n];
	y_data = new double [n];
	//读取数据
	fp=fopen("Nevill_input.txt","r");
	for(int i=0;i<n;i++){
		fscanf(fp, "%lf %lf", &x_data[i], &y_data[i]);
	}
	fclose(fp);//结束读取

	//为了画图，选取x轴范围
	double xmin = x_data[0], xmax = x_data[0];
	for(int i=0;i<n;i++){
		if( xmin > x_data[i] ) xmin = x_data[i];
		if( xmax < x_data[i] ) xmax = x_data[i];
	}
	xmin -= 1;
	xmax += 1;

	int num_grid = 100;//曲线上的点个数
	double step = (xmax - xmin)/num_grid;//曲线上的点的间距

	//计算每个x点的函数值，并输出到文件中
	fp=fopen("Nevill_curve.txt","w");
	if(fp==NULL){
		cout<<"error: failed to open Nevill_curve.txt.\n";
		exit(1);
	}
	for(int i=0;i<num_grid;i++){
		x = xmin + step*i;
		y= Nevill(n,x_data,y_data,x);
		cout<<"y="<<y<<endl;
		fprintf(fp,"%lf   %lf\n", x, y);
	}
	fclose(fp);//结束输出

	return 0;
}
